Introduction to Reservoir Computing with ReservoirPy¶

 Mackey-Glass timeserie¶

Mackey-Glass equation are a set of delayed differential equations describing the temporal behaviour of different physiological signal, for example, the relative quantity of mature blood cells over time.

The equations are defined as:

$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$

where $a = 0.2$, $b = 0.1$, $n = 10$, and the time delay $\tau = 17$. $\tau$ controls the chaotic behaviour of the equations (the higher it is, the more chaotic the timeserie becomes. $\tau=17$ already gives good chaotic results.)

In [2]:
from reservoirpy.datasets.chaos import mackey_glass

timesteps = 25000
tau = 17
X, t = mackey_glass(timesteps, tau=tau)

# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
In [4]:
plot_mackey_glass(X, t, 500, tau)

Prepare the tasks¶

  • Task 1: prediction of the timeserie: from $P(t)$ predict $P(t + 1)$, $P(t + 10)$ (the forecast)
  • Task 2: generation of the timeserie (from $P$ generate $\hat{P}$, and run it over several timesteps).

Task 1 : prepare data¶

In [6]:
forecast = 10  # for now, predict 10 steps ahead
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)

sample = 500
fig = plt.figure(figsize=(15, 5))
plt.plot(X_train[:sample], label="Training data")
plt.plot(y_train[:sample], label="True prediction")
plt.legend()

plt.show()

Task 1: prepare the ESN¶

In [7]:
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
density = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234
In [9]:
im = plt.imread("./static/esn.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()
In [10]:
Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling, 
                                     proba=input_connectivity, input_bias=True,
                                     seed=seed)

W = mat_gen.generate_internal_weights(units, spectral_radius=spectral_radius,
                              proba=density, seed=seed)

reservoir = ESN(leak_rate, W, Win, ridge=regularization)

Task 1: train and test the ESN¶

In [11]:
states = reservoir.train([X_train.reshape(-1, 1)], [y_train.reshape(-1, 1)], verbose=True)

y_pred, states1 = reservoir.run([X_test.reshape(-1, 1)], init_state=states[0][-1], verbose=True)

y_pred = y_pred[0].flatten()
states1 = states1[0]

Computing states: 100%|██████████| 1/1 [00:00<00:00, 2046.00it/s]
Training on 1 inputs (20000 steps)-- wash: 0 steps

Computing states: 100%|██████████| 1/1 [00:00<00:00, 1766.77it/s]
Linear regression...
Linear regression done! (in 0.08027291297912598 sec)
Running on 1 inputs (4990 steps)




















In [13]:



     


sample = 500

fig = plt.figure(figsize=(15, 7))
plt.subplot(211)
plt.plot(np.arange(sample), y_pred[:sample], lw=3, label="ESN prediction")
plt.plot(np.arange(sample), y_test[:sample], linestyle="--", lw=2, label="True value")
plt.plot(np.abs(y_test[:sample] - y_pred[:sample]), label="Absolute deviation")

plt.legend()
plt.show()



     


















    
    



























Determination coefficient $R^2$ and normalized root mean square (NRMSE) :











In [14]:



     


r2_score(y_test, y_pred), nrmse(y_test, y_pred)



     


















    
    
Out[14]:








(0.9999869380130667, 0.0003004561614127586)




















Task 1 : make the task harder¶










In [15]:



     


forecast = 50  # now, predict 50 steps ahead
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)

sample = 500
fig = plt.figure(figsize=(15, 5))
plt.plot(X_train[:sample], label="Training data")
plt.plot(y_train[:sample], label="True prediction")
plt.legend()

plt.show()



     


















    
    



























In [17]:



     


sample = 500

fig = plt.figure(figsize=(15, 7))
plt.subplot(211)
plt.plot(np.arange(sample), y_pred[:sample], lw=3, label="ESN prediction")
plt.plot(np.arange(sample), y_test[:sample], linestyle="--", lw=2, label="True value")
plt.plot(np.abs(y_test[:sample] - y_pred[:sample]), label="Absolute deviation")

plt.legend()
plt.show()



     


















    
    



























Determination coefficient $R^2$ and NRMSE:











In [18]:



     


r2_score(y_test, y_pred), nrmse(y_test, y_pred)



     


















    
    
Out[18]:








(0.9953476933554506, 0.010759732086259792)




















Task 1 : diving into the reservoir¶
Let's have a look at the effect of some of the hyperparameters of the ESN.











Spectral radius¶
The spectral radius is defined as the maximum eigenvalue of the reservoir matrix.











In [19]:



     


reservoir = reset_esn()

states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
    W = mat_gen.generate_internal_weights(units, spectral_radius=sr, 
                                          proba=density, seed=seed)
    reservoir.W = W
    s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
    states.append(s[0].T)



     


















    
    






Computing states: 100%|██████████| 1/1 [00:00<00:00, 2545.09it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 4198.50it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 3045.97it/s]




















In [20]:



     


units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(radii)*100+10+i+1)
    plt.plot(s[:, :units_nb])
    plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()



     


















    
    



























    

  • $-$ spectral radius $\rightarrow$ stable dynamics
    

    • 

  • $+$ spectral radius $\rightarrow$ chaotic dynamics
    

    • 



In most cases, it should have a value around $1.0$ to ensure the echo state property (ESP): the dynamics of the reservoir should not be bound to the initial state chosen, and remains close to chaos.


This value also heavily depends on the input scaling.











Input scaling¶
The input scaling controls how the ESN interact with the inputs. It is a coefficient appliyed to the input matrix $W_{in}$.











In [21]:



     


reservoir = reset_esn()

states = []
scalings = [0.1, 1.0, 10.0]
for iss in scalings:
    Win = mat_gen.generate_input_weights(units, 1, input_scaling=iss,
                                     proba=input_connectivity, input_bias=True,
                                     seed=seed)
    reservoir.Win = Win
    s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
    states.append(s[0].T)



     


















    
    






Computing states: 100%|██████████| 1/1 [00:00<00:00, 807.84it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 4092.00it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 4419.71it/s]




















In [22]:



     


units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(scalings)*100+10+i+1)
    plt.plot(s[:, :units_nb])
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()



     


















    
    



























    

  • $+$ input scaling $\rightarrow$ input-driven activities
    • 

  • $-$ input scaling $\rightarrow$ free activities
    • 



The input scaling can also be used to rescale the inputs and ajust their influences.











Leaking rate¶
The leaking rate ($\alpha$) controls the "memory feedback" of the ESN. The ESN states are indeed computed as:

$$ 
s(t+1) =  \underbrace{\color{red}{(1 - \alpha)} s(t)}_{\text{previous states}} + \underbrace{\color{red}\alpha f(u(t+1), s(t))}_{\text{new states}}
$$
where $s$ is the state, $u$ is the input data, $f$ is the ESN model function, defined as:

$$ f(u, s) = \tanh(W_{in} \cdotp u + W \cdotp s) $$
$\alpha$ must be in $[0; 1]$.











In [23]:



     


reservoir = reset_esn()

states = []
rates = [0.03, 0.3, 0.99]
for lr in rates:
    reservoir.lr = lr
    s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
    states.append(s[0].T)



     


















    
    






Computing states: 100%|██████████| 1/1 [00:00<00:00, 2744.96it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 664.71it/s]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 4306.27it/s]




















In [24]:



     


units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(rates)*100+10+i+1)
    plt.plot(s[:, :units_nb])
    plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()



     


















    
    



























Let's reduce the input influence to see what is happening inside the reservoir (input scaling set to $0.1$):











In [26]:



     


units_nb = 20

plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
    plt.subplot(len(rates)*100+10+i+1); plt.ylabel(f"$lr={rates[i]}$")
    plt.plot(s[:, :units_nb])
plt.xlabel(f"States ({units_nb} neurons)"); plt.show()



     


















    
    



























    

  • $+$ leaking rate $\rightarrow$ low inertia, little memory of previous states
    • 

  • $-$ leaking rate $\rightarrow$ high inertia, big memory of previous states
    • 



The leaking rate can be seen as the inverse of the reservoir's time contant.











 Task 2 : generate predictions¶
During this task, the ESN is trained to make a short forecast of the timeserie (1 timestep ahead).
Then, it will be asked to run on its own outputs, trying to predict its own behaviour.











In [27]:



     


units = 500
leak_rate = 0.15
spectral_radius = 0.48
input_scaling = 1.0
density = 0.1
input_connectivity = 1.0
regularization = 1e-6
seed = 1234

reservoir = reset_esn()

forecast = 1
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)

states = reservoir.train([X_train.reshape(-1, 1)], [y_train.reshape(-1, 1)], verbose=True)



     


















    
    






Computing states: 100%|██████████| 1/1 [00:00<00:00, 2786.91it/s]









    
    






Training on 1 inputs (20000 steps)-- wash: 0 steps










    
    















    
    






Linear regression...
Linear regression done! (in 0.19173073768615723 sec)




















In [29]:



     


start = 0; seed_timesteps = 100; nb_generations = 500

init_inputs = X_test[start:start+seed_timesteps]
X_t = X_test[start+seed_timesteps: start+nb_generations+seed_timesteps]

X_gen, states = reservoir.generate(nb_generations, init_inputs.reshape(-1, 1), return_init=True)
X_gen = X_gen.flatten()

plot_generation(X_gen, X_t, init_inputs, nb_generations, seed_timesteps)



     


















    
    



























 Another example: generate Lorenz butterflies¶
Lorenz attractor is another well-known example of chaotic timeserie generator.


Here, the ESN has to predict 3 covariant values at once: the $x$, $y$ and $z$ coordinates of the Lorenz system.


We can also ask the ESN to predict one value given the two others, and so on. Of course, we can also ask for new generations of the timeserie.











In [30]:



     


from reservoirpy.datasets.chaos import lorenz

timesteps = 25000
X, t = lorenz(timesteps)

# rescale between -1 and 1
X = 2 * (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0)) - 1



     















In [32]:



     


plot_lorenz(X, t, 1000)



     


















    
    



























In [36]:



     


fig = plt.figure(figsize=(20, 20))
plt.subplot(211, projection="3d")
plt.plot(X_gen[:, 0], X_gen[:, 1], X_gen[:, 2], lw=3, label="ESN prediction")
plt.plot(X_t[:, 0], X_t[:, 1], X_t[:, 2], linestyle="--", lw=2, label="True value")
plt.plot([], [], ' ', label=f"$R^2 = {round(r2_score(X_t, X_gen), 4)}$")
plt.plot([], [], ' ', label=f"$NRMSE = {round(nrmse(X_t, X_gen), 4)}$")
plt.legend(); plt.show()



     


















    
    



























Not as easy as Mackey-Glass...











More fun with ReservoirPy and chaotic data: multiscroll attractor, Rabinovich-Fabrikant equations, Hénon map...











In [38]:



     


plot_attractor_examples()



     


















    
    



























 Now the pièce de resistance: real data from the real world¶










In [40]:



     


import requests
import pandas as pd

data = requests.get("https://www.data.gouv.fr/fr/datasets/r/d2671c6c-c0eb-4e12-b69a-8e8f87fc224c").json()
df = pd.DataFrame(data)
df.head()



     


















    
    
Out[40]:











  
    

      
      casConfirmes
      deces
      date
      gueris
      reanimation
      testsRealises
      testsPositifs
      hospitalises
      nouvellesHospitalisations
      nouvellesReanimations
      decesEhpad
      casConfirmesEhpad
      casPossiblesEhpad
    

  
  
    

      0
      191
      3
      2020-03-02
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    

    

      1
      212
      4
      2020-03-03
      12.0
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    

    

      2
      285
      4
      2020-03-04
      NaN
      15.0
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    

    

      3
      423
      7
      2020-03-05
      NaN
      23.0
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    

    

      4
      613
      9
      2020-03-06
      NaN
      39.0
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    

  























French Covid-19 data is available at www.data.gouv.fr.











A new task for the braves : predict the number of hospitalized people¶
We selected 4 interesting features: the number of people hospitalized and new hospitalizations per day, and the number of people in intensive care (réanimation) and new admissions in intensive care per day.











In [43]:



     


plot_dataframe(feat_names, df)



     


















    
    



























In [49]:



     


plot_covid_results()



     


















    
    



























What you could try then:¶
    

  • Playing with the parameters
    • 

  • Some feature engineering and "old school" data science
    • 

  • Hyperopt and ReservoirPy hyperoptimization tools
    • 

  • Add some feedback
    • 

  • Try ensembling methods with several estimators
    • 

  • Contact French government when you are done (and save the world)
    •